Structural and dynamical properties of liquid Si. An orbital-free molecular dynamics 

study 
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Several static and dynamic properties of liquid silicon near melting have been determined from 
an orbital free ab-initio molecular dynamics simulation. The calculated static structure is in good 
agreement with the available X-ray and neutron diffraction data. The dynamical structure shows 
collective density excitations with an associated dispersion relation which closely follows recent 
experimental data. It is found that liquid silicon can not sustain the propagation of shear waves 
which can be related to the power spectrum of the velocity autocorrelation function. Accurate 
estimates have also been obtained for several transport coefficients. The overall picture is that the 
dynamic properties have many characteristics of the simple liquid metals although some conspicuous 
differences have been found. 
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I. INTRODUCTION. 

Molecular dynamics (MD) methods are a useful tech- 
nique to study the properties of liquid systems, and 
the last two decades have witnessed a large increase in 
the application of ab-initio molecular dynamics methods 
(AIMD) based on the density functional theoryi (DFT) 
for the treatment of the electron system. This theory 
allows calculation of the ground state electronic energy 
of a collection of atoms, for given nuclear positions, and 
also yields the forces on the nuclei via the Hellmann- 
Feynman theorem. It allows MD simulations in which the 
nuclear positions evolve according to classical mechan- 
ics with the electronic subsystem following adiabatically. 
Most AIMD methods are based on the Kohn-Sham (KS) 
orbital representation of the DFT (KS-AIMD methods) 
which, at present, poses heavy computational demands 
which limit the size of the systems that can be studied 
as well as the simulation times. However, some of these 
constraints are eased in the so-called orbital-free ab-initio 
molecular dynamics (OF- AIMD) method, which by dis- 
posing of the electronic orbitals of the KS formulation 
provides a simulation method where the number of vari- 
ables describing the electronic state is greatly reduced, 
enabling the study of larger samples (several hundreds of 
particles) and for longer simulation times (tens of ps). 

This paper reports the results of an ab-initio molecu- 
lar dynamics simulation on the static and dynamic prop- 
erties of liquid Silicon (1-Si) at thermodynamic condi- 
tions near its triple point. Si is an interesting mate- 
rial with many peculiar properties which, along with 
its technological importance, has stimulated intensive 
theoretical 2 i 3 i 4 i 5 i 6 i 7 i 8 i 9 and experimental 10 ! 11 ! 12 ! 13 ! 14 : 15 
work. Its high-density forms include the crystalline, 
amorphous and liquid phases with the former two be- 
ing covalently bonded and semiconducting and the latter 
one metallic. At melting Si undergoes a semiconductor- 
metal transition along with a density increase of « 
10% and significant changes in the local atomic struc- 



ture. This evolves from an open one, with a tetra- 
hedral fourfold coordination, to a more compact liquid 
structure with a white-tin-type arrangement preserving 
the local tetrahedrality and with an approximate sixfold 
coordinationiiSii 7 . 

Most theoretical studies on 1-Si have focused on its 
static structural properties and have resorted to classical 
MD simulations in which the liquid system was character- 
ized by effective interatomic potentials constructed either 
empirically by fitting to experimental dataS*^, or derived 
from some approximate theoretical modelA 5 ^. More re- 
cently, some KS-AIMD calculationsZ&Si have been re- 
ported, although the results dealt with electronic and 
static properties only because of the computational con- 
straints inherent in this method. 

Stich et a? performed the first KS-AIMD calculation 
for 1-Si near the triple point using 64 particles, a non-local 
Bachelet-Hamann-Schluter type pseudopotentiali& and 
the local density approximation (LDA) for the electronic 
exchange and correlation energies. Subsequently, a more 
comprehensive study was performed^ with a larger num- 
ber of atoms (350 particles) and a spin dependent gener- 
alized gradient approximation for electron exchange and 
correlation. A KS-AIMD calculation has also been per- 
formed by Chelikowsky et using 64 atoms and a non- 
local, Troullier-Martins type, pseudopotential 19 These 
ab-initio studies have provided accurate descriptions of 
the local liquid structure and valuable insights into the 
the valence electron charge densities, but, by necessity, 
they were too short to address dynamical properties (be- 
sides the self-diffusion coefficient) and these properties 
will now be the main focus of the present report which, 
to our knowledge, is the first ab-initio study on the dy- 
namical properties of 1-Si. 

The static structure factor, S(q), of 1-Si has been 
measured by both neutron scattering (NS)a* and X-ray 
(XR) 11 ! 12 ! 13 diffraction. Although all experimental 5(g) 's 
show a distinctive shoulder at the high-g side of the main 
peak, there are some discrepancies in the positions and 
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heights of the first peak, its shoulder and the second peak. 
Experimental investigations of the dynamic structure of 
1-Si have been hampered by the large adiabatic sound 
speed ( w 4000 m/s), which prevents the use of the in- 
elastic thermal neutron scattering (INS) technique for 
investigating the collective excitations for small g-values. 
Nevertheless, the dynamical structure of 1-Si near the 
triple point has been investigated recently by Hosokawa 
et a ? 14 ' 15 using high-resolution inelastic X-ray scattering 
(IXS) which has no such kinematic restrictions. The mea- 
surements by Hosokawa et al investigated the wavevector 
regions 0.02 A" 1 < q < 1.3 A" 1 14] and 0.02 A" 1 < q < 
2.96 A -1 and several dynamical features already ob- 
served in the liquid alkali metals were found, such as the 
existence of collective excitations up to g-values around 
0.5 q p (where q p is the main peak's position of S(q)), 
along with a positive dispersion in the adiabatic sound 
velocity with respect to the hydrodynamic value. 

Section [H] contains a brief description of the orbital- 
free ab-initio molecular dynamics (OF-AIMD) method, 
and gives some technical details particularly on the elec- 
tronic kinetic energy functional and the local pseudopo- 
tential used to characterize the ion-electron interaction. 
In section lrnl we present and discuss the results of the ab- 
initio simulations for several static and dynamic proper- 
ties which are compared with the available experimental 
data. Finally, conclusions are drawn and possible avenues 
for further study are suggested. 



II. THEORY. 

A liquid simple metal may be treated as a disordered 
array of N bare ions with valence Z, enclosed in a vol- 
ume V, and interacting with N c — NZ valence elec- 
trons through an electron- ion potential v(r). The total 
potential energy of the system can be written, within 
the Born-Oppenheimer approximation, as the sum of 
the direct ion-ion coulombic interaction energy and the 
ground state energy of the electronic system under the 
external potential created by the ions, T4 x t(r, {Ri}) = 

E l =i u (l^-^l) . 

E({Ri}) = £ Z2 +E g [p g (r), V cxt (r, {R t })] , (1) 

where p g {r) is the ground state electronic density and 
Ri are the ion positions. According to DFT, the ground 
state electronic density, p g (r), can be obtained by mini- 
mizing the energy functional 

EW)} = T s [p] + E H [p] + E xc [p] + E cxt [p] (2) 

where the terms represent, respectively, the electron ki- 
netic energy, T s [p] , of a non- interacting system of density 



p(r), the classical electrostatic energy of the electrons 
(Hartree term), 



*»[H//^»-ff. <3) 

the exchange-correlation energy, E xc [p], for which we 
have used the local density approximation and fi- 
nally, the electron- ion interaction energy, E ex t[p], where 
the electron-ion potential has been characterized by a 
first principles local pseudopotential constructed within 
DFTifi. 



E ex t[p] = 




In the KS-AIMD methodji T s [p] is exactly evaluated 
by using a set of roughly N c single particle orbitals, but 
at huge computational expense. In contrast, the OF- 
AIMD approach uses an approximate but explicit density 
functional for T s [p] so that the system is described in 
terms of the single function p(f) replacing the large set of 
orbitalsiiiSS*2i Proposed functionals incorporate the von 
Weizsacker term, 

TwW)} = \J dr\Vp(f)\y P (r), (5) 

plus further terms chosen in order to reproduce correctly 
some exactly known limits. Here, we have used an aver- 
age density model^S, where T s = T\y + T a , 

Ta = ToJ *^(rf /3 - 2a fc(r) 2 (6) 
fc(r) = (2k° F ) 3 J dsk(s)w a (2k a F \r-s\) 

k(r) = (3-7T 2 ) 1 / 3 p(r) a , k° F is the Fermi wavevector for 
mean electron density p e = N e /V, and w a {x) is a weight 
function chosen so that both the linear response theory 
and Thomas-Fermi limits are correctly recovered. Fur- 
ther details are given in Ref. [20|. 

Another key ingredient of the energy functional is the 
the local ionic pseudopotential, v ps (r) , describing the ion- 
electron interaction. It has been constructed from first- 
principles by fitting, within the same T s [p] functional, to 
the displaced electronic density induced by an ion em- 
bedded in a metallic medium as obtained in a KS-DFT 
calculation. Further details are given in Ref. |2(j and 
we merely note that the theoretical framework used in 
this study has provided an accurate description of sev- 
eral static and dynamic properties in the liquid Li, Mg, 
Al, Na-Cs and Li-Na systema 20 ! 22 . 
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III. RESULTS 

OF-AIMD simulations have been performed for 1-Si in 
a thermodynamic state characterized by the temperature 
T= 1740 K and ionic number density pi= 0.05551 A -3 
[Tl| . We have considered 2000 ions in a cubic cell with 
periodic boundary conditions. Given the ionic positions 
at time t, the electron energy functional is minimized 
with respect to p(r) = "0(r) 2 , where V , ( 7= ') is a single ef- 
fective orbital which is expanded in plane waves which are 
truncated at a cutoff energy, Ecut — 12.75 Ryd. The en- 
ergy minimization with respect to the Fourier coefficients 
of ip(r) is performed every ionic time step by using a 
quenching method which results in the ground state elec- 
tronic density and energy. The forces on the ions are ob- 
tained from the electronic ground state via the Hellman- 
Feynman theorem, and the ionic positions and velocities 
are updated by solving Newton's equations, with the Ver- 
let leapfrog algorithm with a timestep of 3.5 x 10 -3 ps. 
In the simulations equilibration lasted 10 ps. and the 
calculation of properties was made averaging over 50 ps. 
For comparison we mention that the previous KS-AIMD 
simulations lasted for 1.2 ps. 0.9 ps. and 1.0 ps. 
, which underscores its limitations for dealing with the 
dynamical properties. 

Several liquid static properties have been evaluated 
(pair distribution function, static structure factor and 
bond angle distribution) as well as various dynamic prop- 
erties, both single-particle ones (velocity autocorrela- 
tion function, mean square displacement) and collective 
ones (intermediate scattering functions, dynamic struc- 
ture factors, longitudinal and transverse currents). The 
calculation of the time correlation functions was per- 
formed by taking time origins every five time steps. Sev- 
eral of these functions also depend on wave vector, which 
for our isotropic system is a dependence on q=| q | only. 

A. Static properties. 

The pair distribution function, g(r), and the static 
structure factor S(q) can both be obtained directly from 
the simulation. The latter is plotted in Fig. ^along with 
the X-ray diffraction data of Waseda et aiH^ and the 
neutron data of Gabathuler and SteebiS. 

The OF-AIMD S(q) has a main peak at q p w 2.59 A" 1 
with a distinctive shoulder at w 3.20 A -1 , the second 
peak is at w 5.65 A -1 and the subsequent oscillations 
are rather weak. As noted earlier there are some dis- 
crepancies among the measured structures in q p , S{q p ), 
and the positions of the shoulder and the second peak. 
The experimental results for q p range from 2.49 A -1 ^| 
to 2.78 A -1 [ljj, the position of the shoulder from 3.15 
A" 1 to 3.40 A" 1 Il0| . and the second peak's posi- 
tion from 5.53 A" 1 [lfto 5.70 A" 1 In all cases 
the smaller estimate is from Takeda's XR data 13 and the 
larger from the NS data of Gabathuler and SteebiS. It is 
worth mentioning that the XR data by Waseda et a<ii*i£ 
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FIG. 1: Static structure factor of 1-Si at 1740 K. Open cir- 
cles: experimental NS data, 3,0 '. Full circles: experimental XR 
data^i*— . Continuous line: OF-AIMD simulations. 

gave values closer to the NS ones, namely q p =2.76 A -1 
with the shoulder at 3.35 A -1 and the second peak's 
position at 5.66 A -1 . Nevertheless, the values from the 
OF-AIMD simulation always lie within the range covered 
by the experimental data. 

The OF-AIMD S(q) accounts for the position and 
height of the main peak, although the amplitudes of 
the following peaks and troughs are somewhat underesti- 
mated. Moreover, the shoulder on the high-g side of the 
main peak, which is the most distinctive feature of the 
measured S(q) for 1-Si, is also reproduced by the simula- 
tion although with a smaller height. A similar shortcom- 
ing is visible in the S(q) from the KS-AIMD simulations 
of Stich et a&£, whereas that of Chelikowsky et gives 
a shoulder of similar height as the main peak. 

Extrapolation of S(q) to q — > and use of the relation 
5(0) = piksTKT gives an estimate for the isothermal 
compressibility kt- A least squares fit of S(q) = so + ^g 2 
to the calculated S(q) for q- values up to 0.8 A -1 yields 
5(0)= 0.04 ± 0.005 and k t = 3.0 ± 0.6 (in 10~ n m 2 
Nw" 1 units) for T = 1740 K. This latter value stands 
remarkably close to the experimental result 23 of 2.8 x 
10~ n m 2 Nw _1 and therefore provides confidence in the 
small-g behaviour of our calculated S(q). 

The simulation directly gives the pair distribution 
function, g(r), whose main peak's position, r p , is usu- 
ally identified with the average nearest neighbor dis- 
tance. The OF-AIMD simulation gives r p w 2.43 A, 
which is slightly smaller than the experimental data of 
2.45 A[Hl|], 2.47 A[l3| and 2.50 AjnJ. An estimate of 
the number of nearest neighbors is obtained by integrat- 
ing the radial distribution function (RDF), 4-7rr 2 p%g{r), 
up to a distance r m taken as the position of its first 
minimumSiS^; the present OF-AIMD study gives r m w 
2.95 A, yielding to a coordination number of « 6.0 atoms, 
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FIG. 2: Bond-angle distribution function, (73(6, r c ). The cut- 
off distance r c is (a) the first minimum of the RDF, r c = 2.95 
A (continuous line) and (b) the covalent cutoff r c = 2.49 A 
(dashed line). 



which reasonably compares with the experimental val- 
ues at T=1735 K ranging from » 6.4 atomsii to w 5.7 
atomsi^. 



Signatures of possible remnants of covalent bonding 
that may remain in 1-Si should be contained in higher 
order correlation functions. A convenient quantity is the 
bond-angle distribution function, g${0, r c ), where 9 is the 
angle between two vectors joining a reference particle 
with two neighboring particles at a distance less than r c . 
In liquid simple metals, this function peaks at around 9 
ps 60° and 120°, which are close to those expected for a 
local icosahedral arrangement (9 « 63.5° and 116.5°). 
The OF-AIMD results for gs(6, r c ), with r c taken as the 
first minimum of the RDF, are plotted in Fig. [21 where 
we see two maxima centered around 9 « 60° and 88°, 
which practically coincide with the values obtained from 
the KS-AIMD simulation*^, namely 9 w 60° and 90°. 
This double-peak feature has usually been interpreted as 
a signature of both tetrahedrally bonded and higher coor- 
dinated atoms contributing to the first coordination shell; 
indeed, according to Stich et afi about 30% of the atoms 
in the first coordination shell are covalently bonded. Also 
plotted in Fig. Elis the g 3 (9,r c ) with r c = 2.49 A which 
has been identified^ as the valence cutoff distance for co- 
valent bonds. Now, the OF-AIMD results show a sharp 
maximum at 9 rs 65° along with broader maximum cen- 
tered around the tetrahedral angle 9 « 109°. Summing 
up, the previous results emphasize the capability of the 
OF-AIMD method to account, albeit not so accurately as 
the KS-AIMD methods, for the orientational correlations 
in those liquid systems with some remnants of covalent 
bonding. 
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FIG. 3: Normalized intermediate scattering functions, F(q,t), 
at several g-values (in A -1 units), for 1-Si at T = 1740 K. (a) 
q= 0.46 (full line), q=0.72 (dashed line), q=1.07 (dotted line) 
and q=1.46 (dot-dashed line), (b) q=2.0 (full line), q=2.57 
(dashed line), q=3.0 (dotted line) and q=3.62 (dot-dash line). 



B. Dynamic properties. 

1. Collective dynamics. 
The intermediate scattering function, F(q,t), defined 

as 

( ? ) 

contains spatial and temporal information on the collec- 
tive dynamics of density fluctuations. Results for F(q, t) 
as functions of time for for several q- values are shown in 
Fig. El F{q,t) oscillates up to?« (3/5) q p = 1.55A" 1 , 
with the amplitude of the oscillations being stronger for 
the smaller g-values. This behaviour, which is plotted in 
Fig. |HIa), is typical of simple liquid metals near melting, 
as found from both computer simulations22i2L22iSi and 
theoretical models^ However, at low q-values (q < 0.5 
q p ) the F(q,t) shows a strong diffusive component which 
plays a dominant role and imposes a slow decay. This is 
at variance with the results for the simple liquid metals 
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FIG. 4: Contributions to the F{q,t) at q =0.47 A" 1 for 1-Si 
at T=1740 K. Diffusive component (dashed line), cosine com- 
ponent (dotted line) and sine component (dash-dotted line). 



near melting)SSiSL22i22^S where for a comparable q-range 
the diffusive component is already very weak and the 
corresponding F(q, t) shows marked oscillations around 
zero. 

This different behaviour can be explained in terms of 
the hydrodynamic expression for the F(q, t), which is ex- 
act in the q — > limit 2 ^& (hydrodynamic limit) 



F(q,t)/S(q) 



1 
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exp(-D T (q)t) 



— exp{—T(q)t) [cos (c s qt) + bq sin (c s qt)} (8) 
7 

and its time FT, the dynamic structure factor, S(q,uj), 



0.015 



3 0.01 
3 

3 0.005 

00 







0.02 



0.01 









(a) 


~\ A 




i 

— 1 l l l 1 l l T— 


\ 

\ 




V 


(b) 


" \ 




\ A 

- \ .a 
















1 1 1 1 1 1 "i • i> 









25 



50 



75 



CO (ps" ) 



FIG. 5: Dynamic structure factor S(k, uj) at several q- values 
(in A" 1 units), for 1-Si at T = 1740 K. (a) q= 0.46 (full line), 
q=0.58 (dashed line), q=0.72 (dotted line) and q=1.07 (dot- 
dashed line), (b) q=1.46 (full line), q=2.00 (dashed line), 
q=2.57 (dotted line) and q=3.00 (dot-dash line). 
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2D T (q) 



7 ) OJ 2 + {D T {q)f 



T(q) + bq{uj + c s q) T(q) + bq(cj - c s q) 



j + c s q) 2 + (T(q)f (w - c s qf + {T{q)f 



(9) 



where 7 = C p /C v is the ratio of specific heats, c s is the 
adiabatic sound velocity, T(q) = Tq 2 with T the sound 
attenuation constant and Dx(q) = Dxq 2 where Z?t = 
kt / (piCp) is the thermal diffusivity and kt is the thermal 
conductivity. Finally, T = ^[a(j ~ 1) /j + vi] with a = 
kt / (PiC v ) and v% is the kinematic longitudinal viscosity. 

According to Eq. JHJl, the hydrodynamic S(q,uj) has 
two (inelastic) propagating peaks centered at u> = ±c s q 
with half- width at half-maximum (HWHM), T(q), and a 
diffusive peak at uj = and whose width is determined 
by the thermal diffusivity, Dt- For a metallic system, 
Dj- has electronic and ionic contributions with the for- 
mer dominating; however it has been shown^i that Eqs. 



©-© incorporate just that part of Dt due to the ions. 
Consequently, we have fitted the OF-AIMD F(q,t)'s at 
the smallest g-values (q < 0.6 A" 1 ), to the analytical 
expression of Eq. JSJ. A good fit was achieved in that 
g-range giving an approximate value Dt ~ 1.0 x 10~ 3 
cm 2 /sec, which is two orders of magnitude smaller than 
the experimental oneS D T = 0.228 ± 0.004 cm 2 /sec. 
for 1-Si at melting, which obviously includes both ionic 
and electronic contributions. We note that an KS-AIMD 
study for liquid Ge near melting 31 gave a comparable es- 
timate of the ionic contribution, namely Dt ~ 1.3 x 10 -3 
cm 2 /sec, which is also two orders of magnitude smaller 
than its total Dt- On the other hand, estimates of the 
ionic contribution to Dt in the liquid alkali metals near 
melting 33 range from 20.0 x 10~ 3 cm 2 /sec. (Li) to 3.0 
x 10~ 3 cm 2 /sec. (Cs). Consequently, as Dt determines 
the diffusive behaviour of the F(q, t) at small q's (see Fig. 
UJ, the larger Dt values of the alkali metals imply a much 
weaker diffusive component which is easily overcome by 
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FIG. 6: Dispersion relation for 1-Si at T =1740 K. Full circles: 
peak positions, ui m {q), from the calculated S(q,ui). Triangles: 
experimental u) m {q), from Hosokawa et al (Ref. 1141) . Open 
squares: peak positions, u)i(q), from the maxima of the calcu- 
lated longitudinal current, Ji(q,u>). Dashed line: Linear dis- 
persion with the hydrodynamic sound velocity, v=3977 m/s. 
Grey (full) diamonds with error lines: experimental (calcu- 
lated) values for the HWHM, V(q), of the inelastic peaks. 



the oscillatory parts of the F(q, t). 

S(q,u)), which is the quantity directly obtained from 
INS or IXS experiments, has been obtained by a time 
FT of the F(q, t) and Fig. [5] shows several S(q, ui) as 
functions of u> for different wavevectors up to « 1.2 q p . 
The S(q,uj) show well defined sidepeaks, indicative of 
collective density excitations, up to q w (3/5) q p which is 
also a common feature of the liquid simple metals 20,26 . 

The dispersion relation of the density fluctuations, 
u> m (q), has been obtained from the positions of the side- 
peaks of S(q,w), and is shown in Fig. HJ1 along with 
uji(q), which is the dispersion relation obtained from the 
maxima of the longitudinal current correlation function, 
Jt(q,v) = u> 2 S(q,uj). Figure El also includes the experi- 
mental oj m (q) data of Hosokawa et ctfi* and a line repre- 
senting the dispersion of the hydrodynamic sound whose 
slope gives the experimental 34 adiabatic sound velocity 
c s = 3977 m/s at T=1753 K. The experimental u) m (q) 
data has a positive dispersion, i.e. an increase of uj m (q) 
with respect to the linear hydrodynamic dispersion rela- 
tion value, with a maximum located around 0.8 A -1 and 
whose magnitude amounts to w 15%. Although there is 
fair agreement between our calculated tu m (q) dispersion 
relation and experiment, the uncertainty at the small q- 
values prevents us from unequivocally declaring positive 
dispersion. A similar positive dispersion has been ex- 
perimentally observed in the liquid alkali metalsSi 3 ^!, 
MgS, AlS and Hg 3 ^. On the other hand, we recall that 
at the hydrodynamic region, the slope of the dispersion 
relation curve gives a g-dependent adiabatic sound veloc- 
ity, c s (q), which in the limit q — > reduces to the bulk 
adiabatic sound velocity. 



Another quantity characterizing the collective density 
excitations is the HWHM, T(q), proportional to the in- 
verse lifetime of the excitations. Figure [S] shows the 
experimental data for T(q), which at low q values is 
oc q 2 (hydrodynamic limit) but departs from it as q 
increases 1 ^*^. The figure also includes the T(q) obtained 
from the simulation, which agrees fairly well with the ex- 
periment. In passing we note that to calculate r(g), the 
F(q, t) obtained from the simulation is fit to an eight- 
parameter analytical expression that interpolates among 
the ideal gas, viscoelastic and hydrodynamic modelsaS. 
This method allows the different contributions to S(q, ui) 
to be disentangled and gives an estimate of T(q) which 
in the present calculations is within an error of « 15 %. 

The transverse current correlation function, J t (q,t), 
is not associated with any measurable quantity and can 
only be determined from computer simulations. It pro- 
vides information on the shear modes and its shape 
evolves from a gaussian, in both q and t, at the free par- 
ticle (q — > oc) limit, to a gaussian in q and exponential 
in t at the hydrodynamic limit (q — ► 0), i.e. 



Jt(q 



0,t) = J- e -9 2 v\t\/m Pi 
(3m 



(10) 



where r\ is the shear viscosity coefficient, (3 — (fegT) -1 
and m is the atomic mass. In both the above limits 
Jt{q,t) is positive, but at intermediate g-values there is 
usually a more complicated behavior with well-defined 
oscillations within a limited q-rangaSiSS^i. J t (q,t) as 
functions of time for 1-Si have been obtained from the 
simulation and are shown in Figure for several q- 
values up to q p . The most striking feature is the ex- 
treme weakness of the oscillations around zero. This is 
in remarkable contrast to the enormous amount of MD 
results for different liquid systems such as hard sphere 
systems 41 ^ Lennard-Jones liquids^^i and simple liquid 
me |- a i 9 20 i 26 J 29 near melting, where the Jt (q, t) exhibit 
strong oscillations and the associated spectra, J t (q,uj), 
have an inelastic peak which appears at low g-values 
and lasts for a finite g-range. Again, 1-Si near melt- 
ing shows a different behaviour with the corresponding 
Jt(q,uj), which is also plotted in Figure Q exhibiting a 
monotonic decreasing behaviour at all g's. Taking into 
account that the inelastic peak in the Jt(q,uj), is asso- 
ciated with propagating shear waves, the results above 
point to their absence in 1-Si near melting. 

An estimate of the shear viscosity coefficient, 77, has 
been made from J t (q,t) as follows jSMSiiii The memory 
function representation of Jt (q, t) : 



Jt(q,z) 



1 

0m 



z H fj(q, z) 



(11) 



where the tilde denotes the Laplace transform, introduces 
a generalized shear viscosity coefficient, fj(q, z). The area 
under the normalized Jt{q, t), gives /3m Jt(q, z = 0) from 
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FIG. 7: (a) Transverse current correlation function J t (q,t) 
at several (/-values (in A -1 units), for 1-Si at T = 1740 K. 
q= 0.46 (full line), q=0.72 (long dashed line), q=1.07 (short- 
dashed line), q=1.46 (dotted line) and q=2.57 (dash-dotted 
line), (b) Same for Jt (q,w). 
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FIG. 8: Self intermediate scattering functions, F s (q,t), at 
several q- values (in A -1 units), for 1-Si at 1740 K. (a) q= 0.42 
(full line), q=0.54 (dashed line), q=0.72 (short dashed line), 
q=0.86 (dotted line), q=1.02 (dot-dash line) and , q=1.29 
(doble dot-dashed line), (b) q= 1.46 (full line), q=2.0 (dashed 
line), q=2.57 (short dashed line), q=3.0 (dotted line), q=3.6 
(dot-dash line) and , q=4.2 (doble dot-dashed line). 



present simulations, this quantity has been obtained us- 
ing 



which values for fj(q, z = 0) can be obtained and these 
may be extrapolated to q = to give the usual shear 
viscosity coefficient, ij. The OF-AIMD simulations give 
7^=0.75 ±0.1 GPa • ps, which fairly compares with the 
available experimental results^ 77 e2 , p =0.58-0.78 GPa • ps. 



2. Single-particle dynamics. 

The most comprehensive information on the single- 
particle properties is provided by the self-intermediate 
scattering function, F s (q,t), which probes the single- 
particle dynamics over different length scales, ranging 
from the hydrodynamic to the free-particle limit. In the 



F s (q,t) = —(^2 e- m ^ t+to) e m ^ to) ) (12) 
3=1 

and Fig. |H1 depicts the results obtained for several q- 
values. It shows the usual monotonic decay with time; 
but comparison with the liquid simple metals (i.e. alkali, 
alkali-earths, Al)22*2£iSi near their triple point shows that 
at similar q/q p values, the F s {q,t) for 1-Si decays much 
faster. As explained later, this fast decay is related to 
the significantly greater self-diffusion coefficient in 1-Si. 

Closely related to F s (q, t) is the velocity autocorrela- 
tion function (VACF) of a tagged ion in the fluid, Z(t), 
which can be obtained as the q — > limit of the first-order 
memory function of the F s (q, t). However, in the present 
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FIG. 9: Normalized OF-AIMD Z(t) for 1-Si at 1740 K (full 
line). The dashed and dotted lines are the longitudinal, Zi(t), 
and transverse, Zt(t), components respectively, as defined in 
Eq. 1151 . and the full circles represent their sum. 
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FIG. 10: Same as in the previous figure but for the corre- 
sponding power spectrum Z{uj). 

simulations it is more easily obtained from its definition 

Z(t) = (v 1 (t)M0))/{vi) (13) 

Figures RjllOl show. respectively, the Z(t) and its power 
spectrum Z(uj) obtained by taking the Fourier transform 
(FT) of the Z(t). Besides the expected diffusive modes 
at low frequencies, the Z(u>) shows a shoulder at to ss 35 
ps -1 , which is close to the optical vibrational frequency 
of covalent Si, and it has been related to the existence of 
vibrational modes echoing the remnant covalent bonds in 
the liquid^. 

The Z(t) shown in Fig. flacks the usual backscatter- 
ing behaviour observed in the liquid simple metals near 
melting which is usually ascribed to the so-called "cage 



effect" due to the shell of nearest neighbors. The rel- 
atively open structure of 1-Si with m 6 nearest neigh- 
bors gives rise to a negligible "cage effect" and Z(t) 
goes monotonically to zero. Indeed, its shape is re- 
markably similar to that from the KS-AIMD calcula- 
tions of Stich et aS^. The self-diffusion coefficient, 
D, is readily obtained from either the time integral of 
Z(t) or from the slope of the mean square displacement 
SR 2 (t) = (\Ri(t) - Ri{0)\ 2 ) of a tagged ion in the fluid, 
as follows 



D = — Z(t)dt; D = lim SR 2 (t)/6t (14) 

pm J t— too 

Both routes lead to practically the same value for 
£>of-aimd = 2.28 A 2 /ps. We know of no experimental 
results for the diffusion coefficients of 1-Si at any ther- 
modynamic state, but note that the KS-AIMD calcula- 
tions of Stich et aU& gave -Dks-aimd = 2.02 A 2 /ps. 
which slightly increased to 2.4 A 2 /ps. when the num- 
ber of particles was increased to 350 particles. They 
found 8 that a spin dependent treatment of electron ex- 
change and correlation within the GGA led to a fur- 
ther increase to -Dks-aimd = 3.1 A 2 /ps., which was 
explained by a weakening of the interatomic bonds in 
comparison with the LDA treatment. The KS-AIMD 
study of Chelikowsky et aU> yielded -Dks-aimd =1-90 
A 2 /ps. The present OF-AIMD result remains within the 
set of values predicted by the KS-AIMD method. Other 
calculations based on tight binding or interatomic pair 
potentials 2 *^ tend to produce significatively smaller val- 
ues for the self-diffusion coefficient along with a VACF 
which takes negative values. We note that ab-initio es- 
timates of the self-diffusion coefficient of 1-Si are nearly 
one order of magnitude greater than corresponding ones 
for the liquid simple metals (alkalis, alkali-earths, Al) 
near meltingSS^^. This larger D explains why F s (q,t) 
decays, for any q-value, much faster than in the simple 
liquid metals. The accurate gaussian approximation^^! 
gives F s (q,t) = exp[—q 2 SR 2 (t)/6], and a greater D im- 
plies a larger 5R 2 (t) and a faster decay for the F s (q,t). 

With knowledge of the self-diffusion coefficient D and 
regarding the motion of an atom as the Brownian mo- 
tion of a macroscopic particle of diameter d, the Stokes- 
Einstein (SE) relation, 77 D = ksT /2nd can be used to 
obtain an estimate for the viscosity rj. Identifying d with 
the position of the main peak of g(r) as is commonly done 
gives d = 2.40 A and taking -Dof-aimd = 2.28 A 2 /ps 
leads to 77=0.69 GPa ■ ps which is consistent with the pre- 
vious OF-AIMD estimate obtained from the transverse 
current correlation function. 

To achieve a deeper physical insight into the features 
of Z(t), we resort to the mode-coupling (MC) approxi- 
mation of Gaskell and Miller— which provides informa- 
tion about the relative importance of the coupling of the 
single-particle motion to the collective longitudinal and 
transverse currents. It has already been used used to 
interpret MD data obtained for the VACF in Lennard- 



Jones fluids, Yukawa fluids, alkali metals and even some 
molecular liquids^^— Within this approach 



Z(t) 



1 



24tt 3 
Z t (t) + Z t (t) 



dqf(q) [Ji(q,t) + 2J t (q,t)}F s (q,t) 



(15) 



where Ji(q, t) and Jt(q, t) are the normalized longitudinal 
and transverse current correlation functions and 



/(«) 



3 ji(aq) 
Pi aq 



(16) 



where j\ (x) is the spherical Bessel function of order one 
and a is an ionic radius given by (4/3)7ra 3 = l/pi- Using 
the OF-AIMD results for Ji(q,t), J t (q,t) and F s (q,t) in 
Eq. 1|15|) allows the contributions of the longitudinal and 
transverse currents to be separated as Zi (t) and Z t (t) re- 
spectively. Firstly we note that, as shown in Fig. the 
Z(t) obtained by applying Eq. (fT5|) is in excellent agree- 
ment with the Z(t) given by the OF-AIMD simulation; 
therefore it may be reliably used to get deeper insight 
into the behavior of the Z(t) by analyzing the separate 
contributions Z[(t) and Z t (t) to the integral (|15|) . 

These contributions to Z(t) are plotted in Fig. Elwhich 
shows that Zi(t) has some oscillatory behaviour, but 
Z t (t) remains positive for all times as a consequence of 
the mostly positive values taken by the Jt (q, t) 's (see Fig. 
0). The positive bump of Z{t) arises from the fast decay 
of Zi (t) which after a negative minimum at a rather short 
time (~ 0.06 ps) quickly goes to zero. This behaviour, 
along with the comparatively slower decay of Z t (t) leads 
to the bump at « 0.12 ps. At longer times, the dynamics 
of Z(t) is completely determined by Z t (t). This shape has 
strong similarities with the MD simulation results for liq- 
uid water at room conditions 49 , although the Z(t) of 1-Si 
always remains positive whereas in liquid water it takes 
negative values after the bump. However, in both sys- 
tems the Z(t) is strongly influenced by Zt{t). In contrast, 
the simple liquid metals show a very different pattern 4 ^ 
with both Zi(t) and Z t (t) oscillating about zero. The 
Z(t) also shows oscillations which, beyond the first min- 
imum, are primarily determined by the Zi(t) which also 
dominates the large t behaviour of Z(t). 

The longitudinal and transverse components of the 
power spectrum are shown in Fig. 1101 The small-w be- 
haviour is completely dominated by the Z^uo), whereas 
the shoulder in the Z{oj) is induced by the maximum in 
Zi(lo) at u> w 40 ps -1 . Note that the absence of a peak 
in Z t (ui) can be traced back to the lack of inelastic peaks 
in J t (g,w) (see Fig. 0), which is in turn an expression of 
the failure of 1-Si to sustain shear waves. 

Again, there are substantial differences between the 
dynamical results for 1-Si and those for liquid simple met- 
als for which both Z t (u) and Zi(lS) have clear maxima 
at a;" 1 and ujY 1 with co r t n < loe < to\ n , where cue is the so- 
called "Einstein frequency" of the metal. 26 Consequently, 



< 



o 

3 

3 

W 2 



. (a) 


\ 

\ 

\ 

\ 

. ■ ■ • • ■ . 

1 1 1 1 1 




// 


: (b) 


// 


// 




// 




// 






i i i i i i i i i 







10 



q(A") 



FIG. 11: (a) Normalized HWHM of S s (q,u), relative to its 
value at the hydrodynamic limit, for 1-Si at T = 1740 K (con- 
tinuous line) and liquid Al near melting (dotted line) . The as- 
teriks are the predictions of the mode-coupling theory and the 
dashed line gives the free-particle limit, (b) Same as above, 
but for the normalized peak value S s (q, OJ = 0), relative to its 
value at the hydrodynamic limit. 



the Z(oj) for liquid metals usually has a maximum at « 
and a shoulder at « which, incidentally, is also the 
pattern exhibited by the MD results for liquid water—. 
In their study of the Z{t) for simple liquids 47 , Gaskell and 
Miller, suggested that the appearance of a peak in Z(u) 
at a frequency well below cue was indirect evidence that 
the liquid system can sustain the propagation of shear 
waves. The present results which show that 1-Si lacks 
both the frequency peak and shear waves, are consistent 
with their suggestion. Finally, we note that the diffusion 
coefficient, which is proportional to Z(lo — 0), is com- 
pletely determined by the transverse component which 
contributes « 99 % of the total. 

The time FT of F s (q,t) gives its frequency spectrum, 
S s (q,uj), known as the self-dynamic structure factor, 
which has experimental importance due to its connection 
with the incoherent part of the INS cross-section. The 
S s (q, to) exhibits, for all q- values, a monotonic decay with 
frequency from a peak value at w = 0. Usually, S s (q,u)) 
is characterized by the peak value, S s (q,u = 0), and the 
HWHM, uji/2(q), which are frequently reported normal- 
ized with respect to the values of the hydrodynamic (q 
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— > ) limit, by introducing the dimensionless quantities 
£(<?) = nq 2 DS s (q,Lj = 0) and A{q) = uj 1/2 {q)/q 2 D. The 
OF-AIMD results for A(q) and E(g) are shown in Fig. El 
along with OF-AIMD results for liquid Al for compari- 
son. Both magnitudes exhibit a behaviour completely 
different from that of liquid Al which, in turn, stands as 
typical of the simple liquids near melting. Indeed, it has 
been founo— that liquid simple metals near their 
triple point have an oscillating A(q), stemming from the 
" cage effect" , whereas for a dense gas it decreases mono- 
tonically from unity at q—0 to a 1/q behaviour at large q. 
The results for 1-Si for both A(q) and T,(q) are more simi- 
lar to those for a dense gas, which must be a consequence 
of the open structure. An indirect check on the reliability 
of these results may be provided by the MC theory£lii£i 
which, has already shown its capability to account for the 
experimental A(q) and S(q) in liquid simple metals^iS 
at q < q p . According to the MC theory 

A(q) = l + H(S)q/q* (17) 
E(g) = l + G(r 1 )< Z /q* 

for small g-values, where q* = \&irmpi(3D 2 , <5 = D/(D + 
•q/mpi) and the functions H(5) and G(8~ 1 ) are given 
in Ref. l52l In both equations the first term is the hy- 
drodynamic result and the second one accounts for the 
coupling of mass diffusion with the collective modes. The 
OF-AIMD results for D and rj have been used to deter- 
mine A(q) and E(g) given by Eq. i|17fl an d the results 
are also shown in Fig. EJ F° r Q < Qp the OF-AIMD and 
MC theory results are in good agreement, similar to that 
found for simple metals. This agreement points to the 
ability of the MC theory to describe the single particle 
dynamics, and presumably the collective dynamics too, 
in 1-Si. 



IV. CONCLUSIONS. 

Several static and dynamic properties of 1-Si at a ther- 
modynamic state close to the triple point have been eval- 
uated. The simulations have been performed using the 
orbital free ab-initio molecular dynamics method com- 
bined with a first-principles local pseudopotential con- 
structed within the same framework. 

The results obtained for the static structure are 
comparable to those obtained by earlier KS-AIMD 
calculations^*^, and are in reasonable agreement with 
the available experimental data. We also stress the good 
description provided by the OF-AIMD method for the 
orientational correlations, as described by the bond-angle 
distribution function, g^,{9,r c ). The obtained structural 
magnitudes are in vein with other previous ab-initio stud- 
ies which suggested for 1-Si a local structure consisting 
in a mixture of both tetrahedrally bonded atoms and a 
higher coordinated structure (probably a distorted metal- 
lic white-tin structure) . This latter assertion is based on 



the idea— that for non-close packed systems, the struc- 
ture of the molten state resembles that of the high pres- 
sure solid state. If we now recall that in crystalline Si the 
semiconducting diamond structure contracts with pres- 
sure and transforms at 12 GPa to the metallic white-tin 
structure, this latter one looks as the most likely to be 
appear in 1-Si. 

Although the ab-initio methods based on the Kohn- 
Sham approach provide a deeper insight into those 
physical magnitudes related to bonding and electronic 
properties^, their heavy computational demands have 
precluded their application to the study of the dynamic 
properties of liquids. Indeed, the calculation of the dy- 
namical properties in 1-Si has been the main aim of the 
present calculations which, incidentally were spurred by 
recent IXS data*^*— Moreover, we stress that this is the 
first ab-initio study on the dynamical properties of 1-Si. 

The intermediate scattering functions, F(q,t), have at 
low-q values, a strong diffusive component which is com- 
parable to what has been found in liquid Ge, but is at 
variance with the liquid simple metals where the diffu- 
sive component is much weaker. The dynamic structure 
factors, S(q, w), show collective density excitations over a 
similar range of wavelenths as those in liquid simple met- 
als. Moreover, the dispersion relation of the excitations 
closely follows the existing experimental data— *— 

The transverse current correlation functions, Jt{q,t), 
show extremely weak oscillations around zero and take 
positive values for most of the time. This shape leads to 
spectra, J t (q,w), with no inelastic peaks which in turn 
reflect the absence of shear waves. This conclusion is rein- 
forced by the analysis of the spectra of the VACF and ap- 
pears as an effect of the low coordination number with its 
attendant negligible "cage effect". This is also reflected 
in the single particle dynamics, as embodied in the VACF 
and the self- intermediate scattering functions, leading to 
a behaviour different from that typical of the liquid sim- 
ple metals near melting. Nevertheless, the MC theory 
appears capable of accounting for the single-particle dy- 
namics in 1-Si. Calculated self-diffusion and shear vis- 
cosity transport coefficients were in fair agreement with 
experiment and/or other ab-initio calculations. 

Up to now, the OF-AIMD method has shown its ability 
to provide an accurate description of the bulk static and 
dynamic properties in simple liquid systems^}*— How- 
ever, the present results for the static and dynamic prop- 
erties of 1-Si illustrate the potential of the OF-AIMD 
method for treating liquid systems with remnants of co- 
valent bonding and therefore open new venues concerning 
its range of applicability. 
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